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The paper presents a numerical study for the finite element method with an- 
isotropic meshes. We compare the accuracy of the numerical solutions on quasi- 
uniform, isotropic, and anisotropic meshes for a test problem which combines sev- 
eral difficulties of a corner singularity, a peak, a boundary layer, and a wavefront. 
Numerical experiment clearly shows the advantage of anisotropic mesh adaptation. 
The conditioning of the resulting linear equation system is addressed as well. In 
particular, it is shown that the conditioning with adaptive anisotropic meshes is 
not as bad as generally assumed. 

1 Introduction 

Anisotropic mesh adaptation, i.e., adaptation of the size and shape of mesh elements, has been 
shown to be of significant advantage for problems with distinct anisotropic features. Moreover, 
the ability to adjust the shape and orientation of mesh elements has proven to be useful for 
designing numerical schemes with particular features, e.g., satisfying the discrete maximum 
principle or improving the conditioning of the finite element equations [71 122] • 

In this paper we concentrate on obtaining anisotropic meshes for the purpose of minimizing 
the numerical solution error. Typically, the optimal shape and orientation of mesh elements 
depend on the Hessian [6l [121 [13 27j or the first derivatives [25l [26j of the exact solution of 
the underlying problem. This is the first major difficulty since the exact solution is usually not 
available. One possibility to solve this difficulty is to try to recover the approximate Hessian 
from the numerical solution in the course of the computation. In \W\ I12|. mesh adaptation is 
based on a residual-based error estimator but still requires Hessian recovery for the solution of 
the dual problem. Unfortunately, Hessian recovery methods do work very well for interpolation 
problems but they cannot provide a convergent recovery if applied to linear finite element 
approximations on non-uniform meshes [3^, '21] , although adaptive finite element methods based 
on Hessian recovery still provide excellent mesh adaptation in practice P [TOl [121 EQl [231 [28] . 

*Supported in part by the National Science Foundation (U.S.A.) through grant DMS-1115118 and the German 

Research Foundation through grants SFB568/3, SPP1276 (MetStroem) and KA 3215/1-1. 
^Department of Mathematics, the University of Kansas |huang@math.ku.edu| 

■'"Department of Mathematics, the University of Kansas (lkamenski@math.ku.edu) 

^Department of Mathematics, Technische Universitat Darmstadt (lang@mathematik.tu-darmstadt.de) 
Center of Smart Interfaces, Technische Universitat Darmstadt 

Graduate School of Computational Engineering, Technische Universitat Darmstadt 



The fact that the convergence of adaptive algorithms employing Hessian recovery cannot be 
proven directly (since the recovered Hessian does not converge to the exact Hessian) explains 
recent interest in anisotropic mesh adaptation based on some kind of a posteriori error estimates 
which do not depend on the exact solution of the underlying problem jH El El SI El [131 El 
[T9] . Moreover, as shown in [191 Sect. 5.3], using error estimates could be of advantage for 
problems exhibiting gradient jumps or similar discontinuities along internal interfaces because 
methods based on recovery of derivatives could result in unnecessarily high mesh density near 
discontinuities. 

For our study we employ the anisotropic mesh adaptation algorithm from [19] which em- 
ploys a globally defined hierarchical basis error estimate (HBEE) for obtaining the directional 
information. In contrast to the recovery-based algorithms, this method adapts the mesh in 
order to directly minimize the a posteriori error estimate and, thus, relies on the accuracy of 
the error estimator but does not require recovery of derivatives of the exact solution. In this 
sense, the algorithm is completely a posteriori. 

Another major concern when using anisotropic meshes is the conditioning of the finite ele- 
ment equations. Generally speaking, an anisotropic mesh is expected to contain elements of 
large aspect ratic[^ and there exists a concern that an anisotropic mesh will lead to extremely 
ill-conditioned linear algebraic systems and this may weaken the accuracy improvements gained 
through anisotropic mesh adaptation. Fortunately, as it has been recently shown in [22j, the 
conditioning of the stiffness matrix with anisotropic meshes is not necessarily as bad as gener- 



ally assumed, especially in 2D. In Sect. 3.2, we will see that even if the condition number of 
the stiffness matrix with an anisotropic mesh is larger than that with an isotropic mesh, the 
accuracy gained through anisotropic mesh adaptation still clearly outbalances the conditioning 
issues, at least for the example considered. 

The outline of this paper is as follows: a brief description of the adaptation algorithm is 
given in Sect. [2] which is followed by the numerical experiment in Sect. [3j The concluding 
remarks are given in Sect. [4j 

2 Discretization and the mesh adaptation algorithm 

We consider a Dirichlet problem for the Poisson equation 

f-An = /, inf2 

1 = 0, on 50 

where O C is a connected bounded polygonal domain. 

For a given triangulation ThoiQ. and the associated linear finite element space Vh C i/o(0), 
the linear finite element solution Uh G V/^ of ([T]) is defined by 

/ Vvh • Vuh dx= fvh dx, Mvh G T^- (2) 
The finite element space Vh and the finite element solution Uh can be written as 

Vh = span{(/)i, • • • ,0n,^J and Uh = ^^l^j^j, (3) 



^In this paper the aspect ratio of a triangular element is defined as the longest edge divided by the shortest 
altitude. For example, an equilateral triangle has an aspect ratio of 2/^/3 ~ 1.15. 
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where is the standard hnear basis function associated with the j-th vertex and riint is the 
number of interior vertices of the triangulation. Substituting ^ into ^ and taking = (pi 
for i = 1, . . . , flint results in the hnear system 

Auh = F, (4) 

where 



Aij = / V(/)j • V(/)i dx and = / 
Jq Jn 



f(j)i dx. 



Note that in order to obtain the finite element solution Uh we need to solve the linear algebraic 
system Q. Thus, the accuracy of Uh depends also on the conditioning of this system which in 
turn is affected by the choice of the mesh. As mentioned in the introduction, there is a concern 
that anisotropic meshes could lead to extremely ill-conditioned linear systems and this may 
weaken the accuracy gained with anisotropic mesh adaptation. In our numerical experiment 
in Sect. 13.21 we will address this issue in detail. 

In order to construct Th (and, thus, the corresponding Vh) we employ the M-uniform mesh 
approach which generates an adaptive mesh as a quasi-uniform one in the metric specified by 
a symmetric and strictly positive definite tensor M = M{x) [17\. The algorithm starts with 
an initial mesh. For every mesh 7^^^^ we compute the finite element solution u^^^^ which is used 
to compute a new adaptive mesh for the next iteration step. The new mesh is generated as an 

(i) (i) 

M- uniform mesh with a metric tensor computed from uj^ . This yields the sequence 

The mesh adaptation process is repeated until the mesh is M-uniform within a given toler- 
ance (see ^19^ Sect. 4.1] for more details). In our computation we use BAMG (bidimensional 
anisotropic mesh generator [16]) to construct anisotropic meshes for a given metric tensor M. 

Typically, the optimal metric tensor depends on the Hessian of the exact solution [HI 
[12} [T7] which is usually unknown. In this study we follow [19j and employ the hierarchical 
basis a posteriori error estimate (HBEE) to obtain the directional information required for the 
metric tensor M^. The brief idea is as follows (see [19] for details). 

If we have an error estimate Zh such that 

\\u - Uh\\ < C\\zh\\- 

for a given norm ||-|| and if it further has the property IlhZh = with being the interpo- 
lation operator associated with Vh (which is fulfilled by the HBEE), than the finite element 
approximation error is bounded by the interpolation error of the error estimate, 

\\u - UhW < C\\zh\\ = C\\zh - IihZh\\. (5) 

Hence, up to a constant, the solution error is bounded by the interpolation error of the error 
estimate and the mesh can be constructed to minimize the interpolation error of z^] the metric 
tensor does not depend on the Hessian of the exact solution. 

In this study, we are concerned with the error measured in the semi-norm, which is 
the energy norm from Q. Therefore, instead of using the metric tensor developed in [19] 
for the error measured in the norm, we construct the metric tensor which minimizes the 
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interpolation error of measured in the semi-norm. In two dimensions the optimal metric 
tensor is given element- wise by 



K 



— \HK{Zh)\ 



' det( I +—\HK{zh)r 
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where Hxizh) denotes the Hessian of the (quadratic) hierarchical basis error estimate on 
element K and is a regularization parameter to ensure that is strictly positive definite. 
ah can also be seen as an adaptation intensity control: uniform mesh has — oo and if 
q:/j ^ the mesh becomes more adaptive. Usually, is chosen so that about half of the mesh 
elements are concentrated in regions where det(M) is large (see [18j for more details on the 
choice of Mk and ah). 

In our computations we employ the globally defined hierarchical basis error estimate since 
it contains more directional information of the solution than localized versions [19l Sect. 5.1]. 
Moreover, it has been shown that local error estimates can be inaccurate on anisotropic meshes 
[8J. To avoid the cost of the exact solution of the global error problem, we use only a few 
sweeps of the symmetric Gauss-Seidel iteration for the resulting linear system until the relative 
difference of the old and the new error approximations is under a given relative tolerance. This 
proves to be adequate for the purpose of mesh adaptation and the computational cost is 
comparable to that of the Hessian recovery: in the tests, the computation of HBEE is about 
twice slower than Hessian recovery. 

Although the validity of the classical hierarchical basis error estimate z^ for the anisotropic 
case is still unclear, theoretical considerations in [15^ Sect. 6.4] and numerical results in [19] 
suggest that the hierarchical basis error estimate is a reliable source of information when a 
mesh is aligned with the solution. 

3 Numerical experiment 

For the numerical experiment we consider a problem in [24j which combines multiple difficulties. 
It is a Dirichlet problem of the Poisson equation 



—/S.U = /, in 

u — on dVt 



(6) 



where is an L-shaped domain — (—1, 1) x ( — 1, 1) \ [0, 1) x (—1, 0]. The functions / and g 
are chosen such that the exact solution u is given by 



u{x, y) = r^/^ sin(26>/3) + tan"^ (^200 ^ (y + 3/4)2 _ 3/4^ 

^ g-1000((x+v^/4)2 + (y+l/4)2) ^ ^-100(y+l)^ 

where r and 9 are the polar coordinates. The solution has 

• a singular gradient at (0, 0) due to a reentrant corner of the L-shaped domain fi, 

• a circular wavefront with the center in (0,-3/4) and the radius of 3/4, 

• a sharp peak at (— \/5/4, —1/4), 

• and a boundary layer along the line y — —1. 

Figure [1] shows the surface and the color plot of a numerical solution. 
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Figure 1: Surface and color plots of the numerical solution. 



3.1 Accuracy of the numerical solution 

First, we compare the accuracy of the numerical solution for Delaunay (quasi-uniform), adap- 
tive isotropic, and adaptive anisotropic meshes. Examples of mesh types are given in Fig. |2j 
We observe that both isotropic and anisotropic adaptive meshes (Figs. [2b| and [2c| respectively) 
have high mesh density in regions with difficulties but the anisotropic mesh (Fig. [2c]) is clearly 
much better aligned with the steep boundary layer and the wavefront. This is the major dif- 
ference between the isotropic and anisotropic adaptation: the isotropic adaptation can provide 
proper mesh density whereas the anisotropic adaptation can provide both proper mesh density 
and proper alignment of the mesh with the anisotropic features of the solution. 

Fig. [3] shows the error of the numerical solution measured in the energy norm |||^x — Uh\\\^ 
which is equal to H^{^}) semi-norm \u — Uh\H^(^Q) for the example considered. The convergence 
plot shows that an anisotropic adaptive mesh requires ca. 200 times fewer elements than a 
quasi- uniform mesh in order to achieve the same accuracy and ca. 10 times fewer elements than 
an isotropic adaptive mesh. In other words, the finite element solution with an anisotropic 
mesh has a 15 times smaller error than an error of the solution on a quasi-uniform mesh with 
the same number of elements and 3 times smaller than the error achieved by means of an 
isotropic adaptive mesh. The asymptotic convergence order of the error in the energy norm is 
the same for all three kinds of meshes: it is 0(A^~^"^)J^ This is expected since we cannot have 
a better convergence order for anisotropic mesh adaptation but can expect a much smaller 
constant when the solution of the problem has anisotropic features. In our test example we 
gain more than one order of magnitude in comparison to quasi-uniform meshes and about one 
half of the order in comparison to the isotropic adaptation. 

Fig. |3] provides also an interesting insight into the behaviour of anisotropic mesh adaptation. 
For very coarse meshes (N < 300) the resolution is not good enough to capture the anisotropy 
of the solution, the mesh is isotropic and has the same error as with isotropic mesh adaptation. 
The interesting part of the plot is between N ^ 300 and N ^ 1000, where the algorithm starts 
to catch the anisotropic features and the error drops quickly. When the anisotropic mesh is 
fine enough to resolve the anisotropy of the solution (N > 1000), the error convergence rate 
reaches the asymptotic state. 



^Note that 0{N °-^) = 0(h) for quasi- uniform meshes in 2D. 
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(a) Delaunay: 2326 elements; max. aspect ratio 2.8. 




(b) Isotropic adaptive: 2321 elements; max. aspect ratio 3.0. 




(c) Anisotropic adaptive: 2316 elements; max. aspect ratio 24.4. 
Figure 2: Mesh examples and 6.6 times close-up views at the reentrant corner. 
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Figure 3: Energy norm of the finite element error vs. number of mesh elements. 



3.2 Condition number of the stiffness matrix 

In this section, we compute the exact condition number (with respect to the ||-||2 matrix norm) 
for the stiffness matrix of the anisotropic finite elements equations and compare it to the 
conditioning of the finite element equations with isotropic adaptive and quasi-uniform meshes. 

The analysis in [22j for the Laplace operator in 2D shows that the condition number can 
be bounded by a term depending mainly on the number N and the largest aspect ratio of the 
mesh elements. In our numerical experiment the maximum aspect ratio is up to 3.8 for quasi- 
uniform and isotropic meshes and up to 37.9 for anisotropic meshes. Thus, the rough estimate 
on the ratio between the condition numbers of the anisotropic and isotropic systems should 
be about 37.9/3.8 ^ 10.0. This is in perfect agreement with our numerical results presented 



in Fig. 4a which show that the condition number of the linear system with anisotropic meshes 
is about one order of magnitude higher than that with the isotropic meshes. Notice also the 
sudden jump in the condition number for the anisotropic case in the range 300 < N < 1000: 
the algorithm starts to catch the anisotropic features of the solution and the maximum aspect 
ratio of the mesh increases quickly as the mesh becomes more and more anisotropic (cp. the 
corresponding error decrease in Fig. [s]). 

Fig. |4a| also shows that that the asymptotic behaviour of the condition number with an- 
isotropic meshes is at most O(A^logA^) which is only slightly larger than 0{N) in the quasi- 
uniform case. The conditioning with isotropic adaptive meshes is also slightly larger than 
0{N) although still smaller than O(A^logA^). Moreover, if a mesh is only locally anisotropic 
(as in our example), a proper diagonal scaling can reduce the conditioning of the stiffness 
matrix so that it is comparable with the condition number in the uniform case (see [22j for 



more details on diagonal scaling). Fig. 4b shows that the asymptotic rate of the conditioning 
of the scaled stiffness matrix is reduced to essentially 0{N)^ which is comparable to that with 
uniform meshes. 
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(b) After diagonal scaling. 
Figure 4: Condition number of the stiffness matrix vs. number of elements. 



4 Conclusion 

Our numerical experiment shows that for problems with anisotropy the anisotropic mesh adap- 
tation is clearly superior to the isotropic one. In our example, at least a half order of magnitude 
could be gained in accuracy by switching from the isotropic mesh adaptation to the anisotropic 
one. The globally defined hierarchical basis error estimate provides good directional informa- 
tion for the anisotropic mesh generation, provided the number of mesh elements is large enough 
to resolve the anisotropy of the solution. It is worth pointing out that the results in Fig. [3] 
present the error of the final numerical solution, i.e., ajter solving the linear system. Thus, 
even if the condition number of the linear system with anisotropic meshes is larger than that 
with isotropic meshes, the accuracy gained through the anisotropic discretization for problems 
with anisotropic features outbalances possible losses due to the numerical accuracy. 
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